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Abstract 

1 The periodic standing wave (PSW) method for the binary inspiral of black holes and neutron 

stars computes exact numerical solutions for periodic standing wave spacetimes and then extracts 
f"^ , approximate solutions of the physical problem, with outgoing waves. The method requires solution 

£S| ■ of a boundary value problem with a mixed (hyperbolic and elliptic) character. We present here a new 

numerical method for such problems, based on three innovations: (i) a coordinate system adapted 
to the geometry of the problem, (ii) an expansion in multipole moments of these coordinates and 
a filtering out of higher moments, and (iii) the replacement of the continuum multipole moments 
with their analogs for a discrete grid. We illustrate the efficiency and accuracy of this method 
with nonlinear scalar model problems. Finally, we take advantage of the ability of this method to 
handle highly nonlinear models to demonstrate that the outgoing approximations extracted from 
the standing wave solutions are highly accurate even in the presence of strong nonlinearities. 
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A. Background 

m : 

The detection and interpretation of gravitational wave signals from inspiralling black holes or neutron stars requires 
a solution of Einstein's equations for the late stages of the inspiral P, 0,0- Much effort is going into the development of 
O" 1 computer codes that will evolve solutions forward in time. For recent progress see [1-0-E1 Such codes will eventually 
Ch . provide the needed answers about the strong field interaction and merger of the binary objects, but many technical 
challenges of such a computation slow the development of the needed codes. This has led us to propose, as a near 
term alternative, the periodic standing wave (PSW) approach. Elements of this approximation have been introduced 
elsewhere 0,0 @j but are most thoroughly presented in a recent paper ^3] that we will hereafter refer to as "Paper 
I." In the PSW approach, a numerical solution is sought to Einstein's equations, not for a spacetime geometry evolved 
from initial data, but rather for sources and fields that rotate rigidly (i.e. , with a helical Killing vector) and that are 
coupled to standing waves. 

Paper I gives the details of how to extract from this solution an approximation to the problem of interest: a slowly 
inspiralling pair of objects coupled to outgoing waves. Paper I also describes the nature of the mathematical problem 
that must be solved numerically: a boundary-value problem with "standing wave boundary conditions" on a large 
sphere surrounding the sources. The differential equations of this boundary value problem are mixed, elliptical in one 
region (inside a "light cylinder") and hyperbolic in another (outside). 

The method of solution in Paper I was straightforward. The differential equations and boundary conditions were 
implemented with a finite difference method (FDM) in a single patch of standard corotating spherical coordinates. The 
equations were solved with Newton-Raphson iteration of a sequence of linear approximations, and a straightforward 
inversion of each linear approximation. The relative simplicity of this approach was useful to demonstrate the basic 
well-posedness and solubility of the problem and to illustrate the important issues of the PSW method, especially 
the "effective linearity" that explains the accuracy of the PSW approximation for the physical solution. The method, 
however, has severe shortcomings. Multipole moments, and hence spherical coordinates are necessary in the wave 
zone for the imposition of outer boundary conditions and for the extraction of outgoing solutions from standing-wave 
solutions. Spherical, and other standard coordinates are, however, not well suited to resolving the relatively small 
sources of the binary. This is especially true if the sources are to be represented by boundary conditions on the outer 
surface of a source, rather than by explicit source terms. The usual technique for handling such problems is coordinate 
patches and interpolation. This would be particularly inconvenient for the PSW computations since standard iterative 
approaches are inapplicable to mixed equations. 
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In this paper we report on an alternative approach, one that has the disadvantage of adding some analytic complexity 
to the problem, and some worrisome features. But it is a method that gives both remarkably efficient results for model 
problems, and a potentially useful new approach to the coupling of moving sources to their radiation field. This new 
method is based on a coordinate system that is adapted to the local structure of the sources and to the large-scale 
structure of the distant waves. Though the PSW computations have been the proximate motivation for introducing 
an adapted coordinate system, the success with this system suggests that its utility may be more broadly applicable. 
Such coordinates, in fact, have already been exploited, even in numerical relativity. "Cadez coordinates" [ill Il2| . 
a carefully adapted coordinate system of this type, was used in much of the work on head-on collisions of black 
holes, and has more recently found to be useful[L| [Tj| for initial data and apparent horizon finding. Like the Cadez 
coordinates, our coordinate systems will reduce to source-centered spherical polar coordinates in the vicinity of the 
sources, and to rotation-centered spherical polar coordinates far from the sources. 

The core of the usefulness of the adapted coordinates is that the field near the sources is well described by a few 
multipoles in these coordinates, primarily the monopolc of the sources, and that the field far from the sources is well 
described only by a few multipoles in these coordinates. A spectral method (that is, a multipole decomposition), 
therefore, requires only a small number of multipoles. We will demonstrate, in fact, that for mildly relativistic sources 
(source velocity = 30% c), excellent results are found when we keep only monopole and quadrupole terms. 

There is, of course, a price to be paid for this. For one thing, there is additional analytic complexity in the set of 
equations. Another difficulty is the unavoidable coordinate singularity that is a feature of coordinates adapted to the 
two different limiting regions. Still, the potential usefulness of the method, and the success reported here have led to 
us treating this approach as the main focus of our computations in the PSW work. 

B. Nonlinear model probelm 

The innovative features of this method present enough new uncertainties that it is important to study this method 
in the context of the simplest problem possible. We use, therefore, the same model problem as in Paper I, a simple 
scalar field theory with an adjustable nonlinearity. We will find it quite useful to set the nonlinearity to zero for 
comparison with the known solution of the linear problem, since many features of our method are unusual even for a 
linear problem. 

For the description of our model problem, we start with Euclidean space coordinatized by the usual spherical 
coordinates r, 9, 0, and we consider sources concentrated near the points r — a, 9 = n/2, in the equatorial plane, and 
moving symmetrically according to <fi = fit and cp = fit + tt. As in Paper I, we seek a solution of the flat-spacetime 
scalar field equation 

f; a ;f3g aP + XF = V 2 * - 9 t 2 * + XF = Source , (1) 

where F depends nonlinearly on 'J. The explicit form of F will be the same as that in Paper I. This will allow 
comparisons with the results of the very different numerical technique in Paper I, and, as in Paper I, allows a very 
useful comparison of the near-source nonlinear solution with an analytic limit. 

We are looking for solutions to Eq. (JJJ with the same helical symmetry as that of the source motions, that is, 
solutions for which the Lie derivative is zero for the Killing vector £ = <9 t 4- Sld^. It is useful to introduce the 
auxiliary coordinate ip = <p — fit. In terms of spacetime coordinates t, r, 9, ip the Killing vector is simply 9 t and the 
symmetry condition becomes the requirement that the scalar field \& is a function only of the variables r, 9 and ip. 
(We are assuming, of course, that the form of the nonlinear term is compatible with the helical symmetry.) It is useful 
to consider the symmetry to be equivalent to the rule 

~nd v (2) 

for scalar functions. In terms of the r, 9, ip variables, Eq. for ^(r, 9, ip) takes the explicit form 

h I m + ^ (-£) + (^r. - *) + w.,..* - w«, ,3, 

that was used in Paper I. 

C. Outline and summary 

The remainder of this paper has the following organization, and makes the following points. In Sec. [H] we introduce 
the concept of adapted coordinates, comoving coordinates that conform to the source geometry near the source and 
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that become spherical comoving coordinates far from the source. A particular system of adapted coordinates, two 
center bispherical coordinates (TCBC), is introduced in this section. Though these coordinates are not optimal for 
computational accuracy, they have the advantage of analytic simplicity and are the only adapted coordinates explicitly 
used in the computations of this paper. Though the TCBC system is relatively simple it is still sufficiently complex 
that that many details of the use of this method are relegated to Appendix [X] 

In Sec. [H] we also discuss the use of these adapted coordinates in a FDM calculation, and explain the computational 
difficulties we encountered in trying to find stable solutions with this approach. These difficulties led us to use a 
spectral type method with the adapted coordintes. In Sec. IIIII we present the fundamental ideas of expanding the 
solution in spherical harmonics of the angular adapted coordinates. In this section we also explain why we are not, 
strictly speaking, using a spectral method since we do the angular differencing by FDM, not by relationships of the 
spherical harmonics. (For background on spectral methods, and an important recent use of spectral methods in 
numerical relativity, see |15|-1 Furthermore, we keep many fewer multipoles than would in principle be justified by the 
number of points in our angular grid. This "multipole filtering" is one of the most interesting and innovative aspects 
of our method. Because the adapted coordinates in some sense handle much of the computation analytically only a 
few multipoles need be kept. In most of the results, in fact, only monopole and quadrupole moments are kept. 

To illustrate a more standard spectral method, we present in Appendix iBla standard spectral approach to the linear 
PSW problem in two spatial dimensions described in TCBC coordinates. The appendix also uses severe multipole 
filtering and serves to demonstrate in a very different, and generally simpler, numerical context the fundamental 
correctness of multipole filtering. 

For the problem in three spatial dimensions, we have found that a special technique must be used for multipole 
expansion and multipole filtering. A straightforward approach would use the continuum multipoles evaluated on 
the angular grid. We explain in Sec. IIIII why this method involves unacceptably large numerical errors, and why 
we introduce a second innovative numerical technique, one that we call the "eigenspectral method." In place of the 
continuum spherical harmonics evaluated on the angular grid, we use eigenvectors of the angular FDM Laplacian. 
These eigenvectors approach the grid-evaluated continuum spherical harmonics as the grid becomes finer but, as we 
explain in this section, the small differences are very important in the multipole expansion/filtering method. Some of 
the details of the eigenspectral method are put into Appendix [UJ in particular the way the FDM angular Laplacian 
can be treated as a self-adjoint operator. 

Section I1VI starts by presenting the details of the model scalar field problems to which we apply the eigenspectral 
method: the choice of the nonlinearity, and the justification for this choice; the manner in which we choose data on 
an inner boundary taken to be the outer surface of a source; the outer radiative boundary conditions; the Newton- 
Raphson procedure for finding solutions to nonlinear problems; and the method by which we extract approximate 
nonlinear outgoing solutions from computed nonlinear standing wave solutions. 

This is followed, in Sec. [VJ by a presentation of numerical results that demonstrate convergence of the method. 
These results show that the numerical methods are quite accurate despite the inclusion of only a very minimal 
number of multipoles. In addition, the power of the numerical method allows us to compute models with much 
stronger nonlinearity than could be handled with the straightforward FDM of Paper I. For these highly nonlinear 
models we confirm the "effective linearity" that was demonstrated in Paper I with less dramatic models: the outgoing 
solution extracted from a standing wave solution is an excellent approximation to the true outgoing solution, even for 
very strong nonlinearity. Conclusions are briefly summarized in Sec. IVII 

Throughout this paper we follow the notation of Paper l[Io|. (A few changes from the notation and choices of 
Paper I are made to correct minor errors of Paper I: (i) the point source delta function is now divided by a Lorentz 7 
factor, as explained following Eq. (|47|l . (ii) The nonlinearity parameter A was used with inconsistent dimensionality 
in Paper I. Here A is consistently treated as a dimensionless parameter, requiring the insertion of a factor 1/a 2 in the 
model nonlinearity of Eq. (|48|l . 

II. ADAPTED COORDINATES 
A. General adapted coordinates 

For the definition of the adapted coordinates it is useful to introduce several Cartesian coordinate systems. We 
shall use the notation x,y,z to denote inertial Cartesian systems related to r,0,4> in the usual way (e.g., z is the 
rotation axis, one of the source points moves as x — a cos (fit), y — asm (Tit), and so forth). We now introduce a 
comoving Cartesian system x, y, z by 

z = rcos8 x — r sin # cos ((f) — fit) y — r sm6 sin ((f) — fit) . (4) 

In this system, as in the inertial x, y, z system, the z axis is the azimuthal axis. We next define the comoving system 
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FIG. 1: Two sets of comoving Cartesian coordinates. 
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in which the azimuthal Z axis is not the rotation axis, but rather is the line through the source points, as shown in 

Fig.n 

Our goal now is to introduce new comoving coordinates x(X, Y, Z), Q(X, Y, Z), Q(X, Y, Z) that are better suited 
to a description of the physical problem, and that allow for more efficient computation. We will assume that the 
coordinate transformation is invertible, except at a finite number of discrete points, so that we may write X, Y, Z, or 
x, y, z as functions of x, O, 

In terms of the comoving Cartesian coordinates, the helical symmetry rule in Eq. ijSJ) takes the form 



V 9y y dx) \ dX dZ 

Our nonlinear scalar field equation of Eq. Q can then be written, for helical symmetry, as 

a 2 * d 2 ^ d 2 ^ 



C^ + XF 



- u 



~ d ~ d 

Z^r - X^ 

dX dZ 



^ + XF = Source . 



dX 2 dY 2 dZ 2 

This field equation can be expressed completely in terms of adapted coordinates in the form 



+ \F = A 
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ee 
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d x de 



2A 



X* 
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2A, 
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a 2 * 



+ S x — + B e —r, - + — + \F = Sources . 
dx <90 <9$ 

It is straightforward to show that the Ay and i?j coefficients here are given by 



A 

^xx 


= v x - 


vx - n 2 A xx 


Aee 


= ve- 


ve - n 2 Aee 


A<i>$ 


= v$- 


V$ - fi 2 A$$ 




= Vx- 


ve - fi 2 A x0 


A x $ 


- Vx- 


v$ - n 2 A x $ 


Ae* 


= ve- 


V$ - ft 2 A e 4> 


5 \ 


= v 2 x 


- n 2 B x 


Be 


= v 2 e 


- n 2 B e 


i?4> 


= v 2 $ 
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(6) 



(7) 



(8) 

(9) 
(10) 

(11) 
(12) 

(13) 
(14) 
(15) 
(16) 
(17) 

Here the gradients, Laplacians and dot products are to be taken treating the X,Y,Z as Cartesian coordinates, so 
that, for example, 

^ - „ dx dQ d X d<d d X dQ 

Vx • ve = — £^ + -4— + — . (is) 

The form of the Ay and £?j terms in Eqs. J2J|-||T7J| are given, for general adapted coordinates, in Eqs. (|A17|I - 1|A27(1 . 
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B. A specific adapted coordinate system: TCBC 




FIG. 2: Geometric basis for the TCBC adapted coordinates. 

Before discussing general features of an adapted coordinate system, it will be useful to give a specific example. 
For that example, we choose a coordinate system x,0,$ that is particularly simple in form, though (as we shall 
discuss below) not the choice that is numerically most efficient. The chosen coordinates are most easily understood 
by starting with the distances n and T2 from the source points, and with the angles 0i, 02 shown in Fig. [21 The 
formal definitions of the adapted coordinates are 



X 

e 



1 ,„ 



(z-a) 



X 2 +Y 2 



Z + a) +X 2 + Y 2 



1/4 



2ZVX 2 + Y 2 



- i#i + 6 2 ) = - tan -1 I — — 

2 ^ ' 2 y Z 2_ a 2_ X 2_ Y 2 / 

tan" 1 (x/Y\ . 



(19) 
(20) 
(21) 



This choice is sometimes called "two-center bipolar coordinates" , hereafter TCBC, and is equivalent to the zero- 
order coordinates used by Cadez [Tll[l^| . 
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FIG. 3: Adapted coordinates in the x, y plane, and three-dimensional coordinate surfaces. 



An attractive feature of this particular choice of adapted coordinates is that the above relationships can be inverted 
in simple closed form to give 



Z = 
X = 
Y = 



2 + X 2 cos 26 + v/(a 4 + 2a 2 x 2 cos 26 + x 4 

1 

2 



-a - X 



1 cos 26 + v/(a 4 + 2a 2 x 2 cos 26 + X 4 ) 



-a 2 - x 2 cos 26 + v/(a 4 + 2a 2 x 2 cos26 + x 4 ) 



cos $ 
sin$ 



(22) 
(23) 
(24) 



The meaning of the \i © coordinates in the x,y plane (the Z 7 X plane) is shown on the left in Fig. |3 a picture of 
three-dimensional Xj @> an d $ surfaces is shown on the right. 
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The geometrical definition inherent in Fig.[2]suggests that the adapted coordinate surfaces have the correct limit far 
from the sources. This is confirmed by the limiting forms Eqs. Ij22 p -I|24 |) for x a. Aside from fractional corrections 
of order a 2 fx 2 the relations are 

Z — > xcos6 X —> x sm © c °s $ Y — > x sm @sin<f> . (25) 

Near the source point at Z = ±a, the limiting forms, aside from fractional corrections of order x 2 /° 2 j are 

2 2 2 

Z -> ±a + cos (29) X -> sin (26) cos $ ? -> ^- sin (26) sin $ . (26) 
2a 2a 2a 

These limits, and Fig. [21 show that near the source point at Z = a the expression x 2 /2a plays the role of radial 
distance, and 26 plays the role of polar coordinate. (Near the source point at Z = —a, the expression x 2 /2a again 
plays the role of radius, but the polar angle is ir — 26.) Notice that both for the near and the far limit, the polar 
angle is defined with respect to the line through the sources, the Z axis, not with respect to the rotational z axis. 

It is clear that our new system has a coordinate singularity at the origin. Indeed, there must be a coordinate 
singularity in any such adapted coordinate system. The switch from the small- x coordinate surfaces, disjoint 2- 
spheres around the sources, to the large- x single 2-sphere cannot avoid a singularity. 

The remaining specification needed is the outer boundary conditions on some large approximately spherical surface 
X = Xmax- For the monopole moment of the field this condition is simply that the field dies off as 1/x- For the 
radiative part of the field we use the usual Sommerfeld outgoing outer boundary condition dtip — —d r ip, approximated 
as dtip — —d x ip- The fractional error introduced by this substitution is of order a 2 jx 2 ■ The Sommerfeld condition 
itself is accurate only up to order (wavelength/7-). Since the wavelength is larger than a, our substitution r — ► % in the 
outer boundary condition introduces negligibly small errors. To apply the helical symmetry we use the replacement 
rule in Eq. © and the outgoing boundary condition becomes 



& [Z—^-X— = )=n r H — + T <P — +T X — , (27) 

d x \ ox dz ) \ oe a® d x > 



where the Ts are given explicitly in Appendix [S] At large x the outgoing condition can be written 

n ( ^ cos6 . - m , n. 2 ,v 

- = fl (cos*- - _ sm*-] (l + 0(aVx 2 )) • (28) 



The correction on the right is higher-order at the outer boundary x = Xmax and can be ignored. The ingoing boundary 
condition follows by changing the sign of the right hand side of Eq. I|27(l or (|28|l . 

The problem of Eqs. © and (|28|l is a well-posed boundary- value problem analagous to that in Paper I^(j- As 
in Paper I, this problem can be numerically implemented using the finite difference method (FDM) of discretizing 
derivatives. The difference between such a computation and that of Paper I is, in principle, only in the coordinate 
dependence of the coefficients (A xx , Aqq,---) that appear in the differential equation and (r e ,---) in the outer 
boundary condition. 



C. Requirements for adapted coordinates 

For the scalar problem, there are obvious advantages of the coordinate system pictured in Fig. EI First, the surfaces 
of constant x approximate the surfaces of constant $ near the sources, where field gradients are largest, and where 
numerical difficulties are therefore expected. Since the variation with respect to 6 and $ is small on these surfaces, 
finite differencing of 6 and <f> derivatives should have small truncation error. The steep gradients in x, furthermore, can 
be dealt with in principle by a reparameterization of x to pack more grid zones near the source points. An additional, 
independent advantage to the way the coordinates are adapted to the source region is that these coordinates are well 
suited for the specification of inner boundary conditions on a constant x surface. Because of these advantages we 
shall reserve the term "adapted" to a coordinate system for which constant x surfaces near the source approximate 
spheres concentric with the source. 

A second feature of the TCBC coordinates that we shall also require in general, is that in the region far from the 
sources, Xj6,<1> asymptotically approach spherical coordinates, the coordinates best suited for describing the radiation 
field. If the approach to spherical coordinates is second-order in a/r, then the outgoing boundary conditions will be 
that in Eq. Jgj). 
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There are practical considerations that also apply to the choice of adapted coordinates. The coefficients of the 
rotational terms in the equation (i.e. , those involving and B t in Eqs. © _ O30 require computing second derivatives 
of the transformation from Cartesian to adapted coordinates. If those relationships are only known numerically these 
second derivatives will tend to be noisy. For that reason, a desirable and perhaps necessary feature of the adapted 
coordinates is that closed form expressions exist for x(2?, j/, z), and Q(x,y,z). (The expression for <£, the azimuthal 
angle about the line through the source points, is trivial.) It is possible in principle, of course, to have the adapted 
coordinates defined without respect to the Cartesians. In the scalar model problem, the coordinates could be defined 
by giving the form of the flat spacetime metric in these coordinates. The nature of the helical Killing symmetry, 
analogous to Eq. © would still have to be specified of course. The choice of adapted coordinates becomes a much 
richer subject in the case of the gauge-fixed general relativity problem that is the ultimate goal of the work; see 



The TCBC coordinates satisfy all the practical requirements of an adapted coordinate system. In particular, the 
functions x(x, j/, z), and 0(x, y, z), as well as their inverses, are all explicitly known in terms of elementary functions. 
Though the TCBC coordinates are therefore convenient, in addition to being well suited to the problem in Eq. 0}, 
they are not optimal. The perfect coordinates would be those for which the constant \ surfaces agree exactly with 
the constant 'J surfaces. This of course is impossible in practice (and, in addition, would not be compatible with 
the requirement that the coordinates go asymptotically to spherical coordinates). We should therefore modify the 
criterion for the "perfect coordinates" to that of having constant on constant \ surfaces for no rotation (Q = 0). 
The TCBC coordinate system, in fact, does satisfy that requirement for the version of the problem of Eq. Q in 
two spatial dimensions with no nonlinearity, as detailed in Appendix Due to this "near perfection" of the TCBC 
coordinates for the linear two-dimensional problem we found that we were able to achieve very good accuracy for that 
case with moderate rates of rotation. 

These considerations suggest that we could achieve an improvement over the TCBC coordinates, by choosing x 
to be proportional to solutions of the nonrotating case of Eq. Q in three spatial dimensions. Since the nonlinear 
case would result in a solution that is known only numerically, we can follow the pattern of the two-dimensional case 
and choose \ simply to be proportional to the solution of the linear nonrotating three-dimensional problem. The 
coordinate that is orthogonal to this x would have to be found numerically, and would therefore be troublesome. But 
there is no need for and x to be orthogonal. We could, therefore, use the TCBC definition of O in Eq. (|20|l . An 
improved set of adapted coordinates, then, would seem to be 



where r^, 0, are the distances and angles shown Fig. El 

In this paper, we shall report only numerical results from the simplest adapted coordinate system to implement, the 
TCBC coordinates. There are two reasons for this. The first is the obvious advantages of working with the simplicity 
of the TCBC case, and the advantage of having simple explicit expressions for all coefficients in Eq. ||SJ). The second 
reason that we do not use the apparently superior adapted coordinate in Eq. (|29[) . is that we do not expect there to 
be an equivalent for the general relativity problem. In that case there will be several different unknown fields to solve 
for, and there is no reason to think that the optimal coordinate system for one of the fields will be the same for the 
others. 



The wave equation in Eq. |(HJ , along with the boundary conditions Eq. (|28Jl , can in principle be solved by imposing a 
X, 0, $ grid and by using FDM. In practice, numerical problems hinder a straightforward finite difference computation. 
Evidence for this is shown in Fig. in which the error (the difference between the computed solution and the analytic 
solution) for the linear outgoing problem is plotted for different locations of the outer boundary x max . As Fig. 0] 
shows, the quality of our solutions was highly sensitive to small changes in grid parameters, such as the location of the 
outer boundary. We attribute these difficulties to the orientation of the finite difference grid at large distance from the 
source. Loosely speaking, the spherical polar grid is "aligned" with the solutions, and errors are distributed evenly 
on the grid. Adapted coordinates become spherical polar at large distances, but the polar axis is aligned with the 
sources, not with the rotation axis. The result may be a nonuniform distribution of errors, which effectively excites 
spurious modes analogous to modes excited inside a resonant cavity. 

An attractive alternative to FDM is to expand $ in a complete set of functions of the angular coordinates and 
$. Since one of our goals is to describe the radiation in the weak wave zone, and since and approach comoving 
spherical coordinates in the weak wave zone, the natural set of basis functions is the spherical harmonics Yi m (Q, $). 



Ref. [13. 
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III. SPECTRAL METHODS WITH ADAPTED COORDINATES 
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FIG. 4: Error in the computed outgoing linear solution as a function of the location of the outer boundary. Results are 
shown both for straightforward FDM in adapted coordinates and for the eigenspectral method, explained in the text, with only 
monopole and quadrupole terms kept. 



In terms of these we would look for a solution of Eq. (JHJ) in the form 



(30) 



even i m 



(The odd is are omitted due to the symmetry of the problem.) 

The possibility of such a spectral method has been introduced in Paper I as a potentially powerful way of dealing 
with radiation from moving sources. The reason for this is that near the source points the field is nearly spherically 
symmetric, and hence can be described with very few multipoles. Far from the source, the contribution from multipoles 
of order t scale as (af2) f , so the radiation field is dominated by the monopole and quadrupole, and again can be 
described with very few multipoles. It is, therefore, plausible that with very few multipoles — perhaps only the 
monopole and quadrupole — the fields everywere can be described with reasonable accuracy. 

In the multipole method, the expansion in Eq. (|30l) is substituted in Eq. © to give 
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(31) 



The next step is to project out ordinary differential equations. This is most naturally done by multiplying by some 
weight function W(x, ©) and by Y^i m i, and by integrating over all and $. The result is our multipole equations 
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(32) 



where Si m is the multipole of the source term, and where 
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The problem with this straightforward approach to multipole decomposition is that the angular integrals needed for 
the projection are very computationally intensive, and the solutions of the differential equations in \ are very sensitive 
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to the values of the a's, /3's, and 7's, that are computed by these projection integrals. These shortcomings do not 
apply to the 2-dimcnsional version of the helically symmetric wave equation. In that case the projection integrals 
involve only a single integration variable, and it proves to be fairly easy to compute accurate angular integrals. 
We present the straightforward 2-dimensional multipole expansion in Appendix [BJ This is meant to illustrate the 
multipole expansion in a particularly simple context, but more important it demonstrates a crucial point, that we can 
get excellent accuracy by keeping only two multipolcs. This 2-dimcnsional computation also illustrates the alternative 
definiton of standing waves, that of minimum wave amplitude, as sketched in Paper I. 

It turns out that for the 3-dimensional problem, even with only a small number of multipoles, there are two classes 
of severe computational difficulties. First, the projection integrals in Eq. il.'SMl are very computationally intensive, 
especially due to the singularity at 6 = 7r/2 for x/a < 1, a singularity that must be canceled in the projection 
integrals by the choice of the weight function W(x>®)- In trials with the linear problem, and in comparisons with 
the known exact answer, we have found that accuracy of the computed field is poor unless the integrals are done very 
precisely. A second, quite distinct, difficulty is related to the projection at the outer boundary. An outgoing boundary 
condition is applied to d£ m , for I > 0. The radiative moments, however, are much smaller than the monopole moment 
aoQ. Projection of a ae m moment with I > will be contaminated by the much larger monopole moment aoo, due 
to small numerical inaccuracies in the projection. We have found this to be a problem even in the simplest (static 
linear) models. 

We have used an alternative approach to multipole decomposition and multipole filtering, an approach that gives 
excellent results for the nonlinear scalar models and promises to be similarly useful in gravitational models. Underlying 
this approach is the concept that the angular nature of the multipole components of the radiation field is determined 
by FDM operations, in particular by the FDM implementation of the Laplacian. The properties of the spherical 
harmonics that make them useful in the continuum description of radiation is taken over, in FDM computations, by 
the eigenvectors of the FDM Laplacian. 

To implement this idea we start by viewing the grid values of the scalar field $ on a constant-x surface as a vector 
\& whose components are most conveniently expressed with a double index 

*<:<-),. <1\!. (34) 

Here 0; and are the values on the O, <& grid with spacings A0 and A$. It follows that ^ is a vector in a space 
of dimension N = uq x n$. 

In the 0, $ continuum, the angular part of the Laplacian at \j a 3> 1 is the operator 



V 2 EE 

sin de 



ana 



smQ 00 



1 d 2 

(35) 



(sin©) 2 <9$ 2 ' 



In a FDM this is replaced by an operator in the TV-dimensional space of angular grid values. Our eigenspectral method 
is based on finding the eigenvectors of this iV-dimensional operator. 

(k) 

The i,j component of the eigenvector will have the form which should be a good approximation to some 
Yim(®i,&j), i- e - 1 to some continuum spherical harmonic evaluated at grid points. (In practice we work only with 
real eigenvectors that are approximations to normalized real and imaginary parts of the grid-evaluated spherical 
harmonics.) In Fig. [3] continuum spherical harmonics are compared to the eigenvectors found for a grid with ne x 
n$ = 16 x 32 on an angular domain < < n/2, < 4> < ir. As might be expected, the agreement between eigenvector 
and continuum function is quite good when the scale for change of the continuum function is long compared to the 
spacing betwen grid points. 

The eigenvalues found for the discrete and continuum angular Laplacians are in good agreement for small eigen- 
values. For the discrete problem we define an effective multipole index £ in the obvious way, by setting —£(£ + 1) 
equal to the eigenvalue for each eigenvector. A comparison is given in Fig. [5] of the integer continuum values of I 
and those found for a 16 x 32 grid on the region < < tt/2, < $ < ir. (Unlike the spherical harmonics, the 
eigenvectors are not degenerate, so there is a small range of I values of the eigenspectral method corresponding to 
each I of the continuum problem.) For the discrete operator the eigenvectors) For our problem the other angular 
regions are related by symmetry. These symmetries also eliminate the odd values of I omitted from Fig.[|J] The figure 
shows that for small I there is good agreement between the discrete and continuum eigenvalues. Because of this we 
can refer to monopole, quadrupole, hexadecapole, . . . eigenvectors without ambiguity. 

In the mathematics of the grid space, two vectors i*y and Gij are taken to have an inner product 

F ' G = J2 £ F « G « sin(0i)A0 A$ . (36) 

i=l 3 = 1 
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FIG. 5: The eigenvectors for a 16 x 32 grid compared to the corresponding continuum eigenfunctions, the spherical harmonics. 
The continuous curves show the spherical harmonics; the data points are the components of the eigenvectors. 



With some care, detailed in Appcndix|0 we can construct the FDM angular Laplacian to be self adjoint with respect 
to the inner product in Eq. I|36[l . This guarantees that the eigenvectors can be chosen to be orthogonal. We complete 
the analogy to the spherical harmonics by choosing the eigenvectors to be normalized, so that we have 



(37) 



EE^ fe) 4 fc ' )sin (e 4 )AeA<i> = 4 fc '. 
i=l j=l 

(This normalization has been used for the eigenvectors shown in Fig. 0) With these definitions we can now write a 
multipole expansion as 



*(x,e i ,$ J -) = E a(fe) W^ 



(38) 



with 



(x)=EE*(*> < 



,^)Y^ sin(eOAeA$, 



(39) 



= 1 3=1 



The multipole filtering that was the motivation for the introduction of the spectral decomposition is implemented 
simply by limiting the terms included in the sum in Eq. 1)38(1. Rather than include all eigenvectors, only those with 
i < ^max are included. Since the discrete is are never larger than the continuum is, a choice € max — 5 means that the 
monopole, quadrupole and octupole terms, with I ss 0,2,4 are included. The effect of the eigenspectral method and 
multipole filtering are the suppression of the large FDM boundary-related errors, as illustrated in Fig. 21 

In this method the k summation in Eq. (|38|l stops at some maximum value governed by £ max , and the equations to 
be solved are the following modifications of Eqs. I|32|) : 



u d X 2 



Sk' 



In place of Eq. (|33|l the coefficients in this sum are now evaluated from 



OLk'k 



k'k 



Ik'k 



= y( fe/ ) 



A 



ee ■ 



2A 



A 



xe- 



80 



2A 



Q2y{k) 

gy( k ) 



2A« 



e<i> 



Q2 Y (k) 

d¥dO 



gyi k ) gy( fc ) 
-Be — ~ — I - B$ 



d® 



x*" 



9$ 



B x Ye, 



(40) 

(41) 
(42) 

(43) 
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'continuum 



FIG. 6: The i values of the discrete angular Laplacian on a 16 x 32 grid compared with the integer I values of the continuum 
angular Laplacian. The eigenvectors of the discrete angular Laplacian are not degenerate, so a cluster of several I values of the 
eigenspectral method corresponds to a single t value of the continuum problem. 

where it is understood that the angular derivatives are computed by finite differencing. In the effective source term, 

S y = -XY^ ■ F (J^ aWy( fc )) , (44) 

only the nonlinearity appears. There is no "true" source term since we solve only outside the source and introduce 
the properties of the source through boundary conditions. 

Our method is clearly spectral in flavor, but it is worth pointing out explicitly that this method is not a spectral 
method according to the meaning usually given to that term in numerical analysis. If it were a spectral, or pseu- 
dospectral (collocation) method, then angular derivatives in the field equations and boundary conditions would be 
taken using properties of the spectral functions. (If the decomposition were done into continuum spherical harmonics, 
for example, a spectral method would evaluate d^f/dQ by using relations among the spherical harmonics.) In our 
method, angular derivatives are taken by finite differencing, not by relations among the eigenvectors and their angular 
derivatives. We could, in principle, convert our method to one that meets the "spectral method" (actually pseudospec- 
tral) definition. We could use finite differencing to compute, once and for all, relations among the eigenvectors and 
their derivatives. These relations could then be used to replace any derivative by a linear combination of eigenvectors. 
We have, however, not explored this approach. 

Some comment must be made about a subtle but fundamental point in our spectral method. For a given x < a, the 
angular specification © = ir/2 refers to a single point on the Z axis; the value of <!> is irrelevant. On the other hand, 
the function Y^ m (7r/2, $), for even i and is, in general, not a single value. There are, then, terms in Eq. i|3U|) 

that in principle are multivalued at \ < a >© = ir/2. We can, of course, delete the value O = ir/2 from our grid. (And 
we, in fact, delete this value for several reasons, such as the requirement that the FDM Laplacian be self-adjoint; see 
Appendix El) We still have the problem that the variation of these awkward terms diverges as AO — > and the grid 
converges to the continuum. In principle, for any AO the summation in Eq. (|3U[) at any grid point will approach (in 
the mean) the correct answer if we include enough multipoles. 

In practice, we include very few multipoles. We must therefore ask whether the summation will give a highly 
inaccurate answer in the region of the x < a grid near O = tt/2. We avoid this problem by choosing source structures 
that are symmetric about the Z axis. This means that at some inner boundary Xmin we set the nonaxisymmetric 
to zero. The radial equations, the FDM eigenspectral versions of Eq. H32fl . do mix the a^ k \ so the nonaxisymmetric 
a( fc ) will be generated. But the mixing of the multipoles is small until \ is on the order of a. As a consequence, the 
nonaxisymmetric a^ k ' can play their needed role in the wave region without generating large errors in the near-source 
region. 

This behavior of the coefficients is illustrated in Fig. \7\ for an outgoing linear wave. The solid curve shows, as a 
function of x, the eigenspectral coefficient corresponding to the 1 = 2 mode that is symmetric about the source 



12 



axis that is, the mode corresponding to Y2o(0,<&); the dashed curve shows the eigenspectral mode corresponding 
to the real part of Y 2 2(Q,^)- In both cases, the value of the coefficient is divided by the value of the monopole 
coefficient to give a better idea of the relative importance of the mode in determining the overall angular behavior. 
The Y2o(6,$) mode, which does not involve multivalued behavior on the Z axis has a nonnegligible coefficient at 
small x- By contrast the 122(0,$), which is multivalued, has a very small coefficient, one that is two orders of 
magnitude smaller than the monopole, up to x ~ 1- F° r larger x this mode gets "turned on," as it must, since it is 
part of the radiation. 




FIG. 7: The \ dependence of the eigenspectral mode coefficients. The solid curve shows the coefficient of the mode Y20 that is 
symmetric about the Z axis; the dashed curve shows the real part of Y22- In both cases the plot shows the coefficients divided 
by the monopole coefficient. 



IV. MODELS AND METHODS 



Nonlinear scalar models 



The model problem of Paper I, in the original comoving spherical coordinate system is 

£(¥)=0-eff[«l, 

with L taken to be 



1 3 / J 



1 d ( . B 



1 



r 2 dr \ dr J r 2 sin 9 89 \ 89 
and with the effective source terms 

CT c ff [vP] = point source — XF 



r 2 sin 2 1 



n 2 



d 2 

dip 2 



(45) 



(46) 



(47) 



In Paper I an explicit delta function term was used in a e g to represent the point source. Here we compute only outside 
the source and include source effects by the inner boundary conditions described below. 
Our choice of the nonlinearity function F is 



F 



1 * 5 



(48) 
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in which Wo is an adjustable parameter that we set to 0.15 or 0.01 in the numerical results to be reported. As detailed 
in Paper I, this choice of F allows us to make useful estimates of the action of the nonlinearity. We briefly review this 
feature here. 

We define R to mean distance from a source point, and we identify Rn n as the characteristic distance separating the 
|w| ^> | Wo | near-source nonlinear region, and the |w| <C |Wo| distant region in which nonlinear effects are negligible. 
In the nonlinear region near a source of strength Q/a the solution approximately has the Yukawa form 

W sa — — near source pt . (49) 

a 4nR/a 

We can estimate Ru n by taking it to be the value of R at which the expression in Eq. I|49|) is equal to Wo: 



a AirRnn/a 



*o • (50) 



If Rn n is significantly less than a, which it is for most of the models we consider, then we can approximate W as having 
the Yukawa form in Eq. (|49|l out to Ru n . For R > Ru n the linear Coulombic form should apply. We can therefore 
view exp (— \J — A R\ in /a) as a factor by which the strength of the source is reduced. Since the waves are generated at 
distances from the source much greater than Ry m , the wave amplitude as well as the monopole moment of the source 
should be reduced by this factor. We saw in Paper I that these estimates were in reasonably good agreement with 
the results of computation, good enough to give confidence of the fundamental correctness of the picture on which 
the estimate is based. We therefore use this picture in the present paper in interpreting some of the computational 
results. 



Boundary conditions 



In Paper I the source was taken to be two unit point charges moving at radius a 



S = 7 " 1 ^- r ^ 6(0 - tt/2) [%) + <% - tt)] . (51) 



Here 7 is the Lorentz factor 1/yl — a 2 fl 2 . This factor is necessary if the source is to correspond to points of unit 
strength as measured in a frame comoving with the source points. (This factor was inadvertently omitted from the 
source in Paper I. In that paper only the case afl = 0.3 was studied, so we may consider the point sources in Paper I 
not to have been unit scalar charges, but source points with charges 7 = 1.048.) In the present paper we specify inner 
boundary conditions on some surface Xmin rather than an explicit source term as in Paper I. Our standard choice for 
the inner boundary conditions will be those that correspond to the point sources of Eq. (|51|l . For this choice of source 
and for Xmin *C a and <C (— A) _1 / 4 a, we can use an approximation for a single source point. 
In notation appropriate to the 3D case we have 

R 2 = Z 2 +-/ 2 (X - /3t) 2 +Y 2 (52) 

Now we use the transformations of Eq. (|22[1 - l|24|) to get 

R 2 = + ( 7 2 - 1) (x - (it) + 0( X 6 ) = \ [1 + (7 2 - 1) sin 2 28 cos 2 $] (53) 
For a unit strength source at position 1, the field near position 1, due to source 1, should be 

* = H"r4 ; <»«> 

The outer boundary condition in our computation is based on the the Sommerfeld condition in Eq. (|28() ; the ingoing 
condition is identical except for a change of sign. These radiative boundary conditions should be applied only to the 
radiative part of the wave. This is done by applying the conditions to the sum on the right side of Eq. 138|) with the 
monopole mode omitted. The multipole components of this outer boundary condtion are then projected out. The 
monopole moment, of course, is nonradiative. Since it falls off at large distances as 1/x, the outer boundary condition 
is taken to be 



\ d\ X 



= 0. (55) 

Xmax 
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Extraction of an outgoing approximation 



In Paper I, we explained how to extract a good approximation of the outgoing solution from the computed standing- 
wave solution. That explanation started with the solution of the linearized problem 

*stndcomp = 51 ^ a im(r)Y im (6,(p) . (56) 

even^ m=o,±2,±4.. 

We keep that notation here, but understand that (i) the role of the continuum spherical harmonics is played by 
appropriate linear combinations the eigenvectors, that (ii) the role of the coefficients ai m is played by appropriate 
linear combinations of the coefficients a^ k \x): an d that (hi) the summations only extend up to £ max . 

As in Paper I, this form of the computed standing-wave solution is compared with a general homogeneous linear 
(A = 0) standing-wave (equal magnitude in- and outgoing waves) solution of, with the symmetry of two equal and 
opposite sources: 

*tndlln = Y lm (9,cp) [±C im h£\mnr) + ±C* m h¥\mnr)\ . (57) 

even £, m 

A fitting, in the weak-field zone, of this form of the standing- wave multipole to the computed function a^ TO (r) gives 
the value of Ci m . 

By viewing the linear solution as half-ingoing and half-outgoing we define the extracted outgoing solution to be 

*exout = J2 Y tm {9,i P )C em h' e 1) (mnr) . (58) 

eveni m=o,±2,±4.. 

Since this extracted solution was fitted to the computed solution assuming only that linearity applied, it will be a 
good approximation except in the strong-field region. In the problems of interest, the strong-fields should be confined 
to a region near the sources. In those regions, small compared to a wavelength, the field will essentially be that of 
a static source, and will be insensitive to the distant radiative boundary conditions. As explained in Paper I, the 
solutions in this region will be essentially the same for the ingoing, outgoing, and standing-wave problem. In this 
inner region then, we take our extracted outgoing solution simply to be the computed standing-wave solution, so that 

T _ J E^mftm''/ weak field outer region ftiQ \ 
^exout — \ ~r , c i , • • • l oy J 

I ^stndcomp strong field inner region 

The transition between a strong field inner region and weak field outer region can be considered to occur in some 
range of \. The maximum \ i n this range must be small compared to the wavelength l/f2, and the minimum x must 
correspond to a distance from the source larger than our estimate of i?n n - [See Eq. H50|) ]. For distances R from the 
source that are of order a or less, x ~ v2ai? so the the minimum x. in the transition region should be larger than 
X ~ \/2ai?iin ■ 

In order for the extracted solution to be smooth at this boundary, we construct our extracted solution by using a 
blending of the strong-field inner solution and the weak- field outer solution over a range from xiow to Xhigh- In this 
range we take 

^exout = /%) Y trnC'e m h ( p + [1 - f3( X )] *st„dcomp • (60) 



Here 

X-Xlo 



0{X) 



Xhigh Xlow 



X-Xlo 



Xhigh Xlow 



(61) 



so that j3(x) g° es from at x = Xlow to unity at x = Xhigh and has a vanishing x-derivative at both ends. In principle 
we should choose x = Xhigh to depend on the location of the wave zone, and hence on tt, and in principle we should 
choose x = Xiow to depend on the nature of the nonlinearity, and hence on A and "to- In practice we have found it to 
be adequate to choose Xhigh = 2a and Xhigh = 3a for all models. 
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Nonlinear iteration 



The computational problem of finding a solution 'J consists of finding a set of coefficients (x) that satisfy the 
field equation Eqs. I|4(J|I along with the inner and outer boundary conditions. The operations on the left hand side 
of Eqs. I|4l)|l are linear on the aW(x), as are the boundary conditions, so the problem of finding the a^ k \x) can be 
written as 

EW k) =^({« W }), (62) 
fc 

where Ck'k is a linear differential operator on the a^(x)j an d where Tk'({a^}, containing the nonlinearity in the 
model, is nonlinear in the a^(x)- 

For different boundary conditions (outgoing or ingoing) the linear operator Ck'k has different forms, but in either 
form we can invert to get the outgoing or ingoing Green functions >C^,] c ' out and C^, 1 ^ 111 . In principle we can then find 
solutions by direct iteration 

a£lout = E^ out (*v({«SL») C\m = E^ in (^({«SSL») 

k k 



,(*') 



E I {^ out +^ out } (^({«S tnd })) ■ 



(63) 



In Paper I it was pointed out that this kind of direct iteration converges only for weak nonlinearity. More generally 
we use Newton-Raphson iteration and solve 



E 



c^k — 



Ak') 



\ = r» (K p) })-E 



dTy 



l n+l 



(64) 



This Newton-Raphson approach can be applied to find outgoing, ingoing and standing-wave solutions analogous to 
those in Eqs. I|63[) . It has been applied with an error measure 



\ h — ^ ( a n+ite) - 4+1 (Xi) 



(65) 



Iteration was halted when this error measure fell below ~ 10~ 6 . Note that for strongly nonlinear models, convergence 
sometimes required that the iteration described in Eq. (|64|) . had to be somewhat modified. The last term on the right 
in Eq. (|64p. had to be weighted by a factor less than unity, at least until the iteration got close to the true solution. 



V. NUMERICAL RESULTS 



If numerical results are to be trusted they must converge, or at least be stable, as computational parameters (grid 
size, etc.) change, and there must be evidence that the result is the correct answer to the physical problem. A 
complication in demonstrating this is that at the same time we are making two different classes of approximations: 
(i) we use values on a grid in place of the continuum mathematics, (ii) we are keeping only low order multipoles. 
In addition, to represent point sources we use approximations for inner boundary data that are exact only only for 
Xmia ~ * 0- Our outer boundary conditions in Eq. I|55f) also add an error, in principle one of order (a/xmax) 2 , but we 
have found that this error is negligible compared to that of our other approximations. (Moving the boudary outward 
has no discernible effect on results.) Here we present results of varying the grid resolution, the number of multipoles 
kept, and the inner surface Xmin on which inner Dirichlet data is set. 

Error is most easily measured for solutions to the linear problem since there exists an exact series solution for 
comparison. Figure [S] shows a comparison of this series solution with computed solutions for two different source 
speeds afi. The qualitative features of these plots agree with what should be expected: the eigenspectral/multipolc 
filtering technique is more accurate at lower source speed, and is more accurate when more multipoles are allowed 
to pass through the "filter." (Of course, there will be a point of diminishing returns. If we let too many modes 
through then we are no longer filtering, and we experience the difficulties that plagued the FDM method for adapted 
coordinates.) 
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r/a r/a 



FIG. 8: Comparison of exact and eigenspectral linear outgoing solutions. The solid curve shows the exact solution, in the 
wave region, computed from an infinite series. The other curves show the result of computation with a grid with n x x n& x 
n$ = 12001 x 16 x 32 and Xmin/o = 0.2, Xmax/a = 75. Results are shown with l m&K = 3 (monopole and quadrupole modes kept) 
and i'max = 5 (monopole, quadrupole and hexadecapole) . The results are shown at O = as a function of r, the distance from 
the center of the configuration. 

TABLE I: Convergence for rotating linear models. All models have ail = 0.3, A = 0, and use outgoing boundary conditions 
at Xmax = 50a. The computed monopole to source strength index, jQad/Q, is unity in the exact solution. The "two region" 
computation retains all multipoles for \ < 3a. 



n x 


ne 


71$ 


Xmin/a 


p 
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yQea/Q 


two region 
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1.0116 
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1.0275 
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1.0282 
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1.0284 
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1.0282 
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0.2 
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1.0036 


1.0035 


8001 


16 


32 


0.2 
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0.9968 


0.9966 


8001 


16 


32 


0.2 
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0.9934 


0.9928 


16001 


16 


32 


0.2 


5 


1.0037 




16001 


16 


32 


0.1 
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1.0036 




16001 


16 


32 


0.05 


5 


1.0032 




16001 


16 


32 


0.025 


5 


1.0007 





An obvious unwelcome feature of the results with £ max = 5 is the phase error in the hexadecapole mode. This 
error is especially noticeable in the a£l = 0.5 plot where it causes an artifact fine structure at the positive peak of 
the waves. We are investigating the source of this phase error which we suspect is a result of truncation error in 
angular differencing and/or in the computation of the eigenvectors. We have anecdotal evidence that the phase error 
decreases as the angular grid is refined. 

The reliability of the eigenspectral method for a wider variety of linear models is presented in Tabled I n this table, 
the measure of error is the value of Q e ff, the monopole moment computed near the outer boundary for the "charge," 
i.e. , the monopole moment of the two sources each with scalar charge Q/a = 1. Though the monopole moment would 
seem to be less interesting than features of the radiation, we have found in essentially all computations that the largest 
error is in the monopole. For example, the majority of the error in the computed amplitude of radiation could be 
understood to be due to the error in the monopole. The error in the ratio of radiation amplitude to monopole was 
several times smaller than the raw errors in either quantity by itself. For simplicity we use this one measurement to 
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TABLE II: Convergence for rotating nonlinear models. All models have aQ = 0.3, A = — 25, 'I'o = 0.15, and all use outgoing 
boundary conditions at Xmax = 50a. The "two region" method retains all multipoles for \ < 3a. 
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characterize convergence and correctness. 

To show that the computed solution is accurate it is convenient to consider first the linear outgoing problem for 
two unit point charges, since the solution for this case is known to be Q c s — I/7, where 7, the Lorentz factor, is 
1/Vl — a 2 il 2 . (The complete solution for in this case is given as Eq. (10) of Paper I, though the series solution 
must be multiplied by I/7 since we are now considering unit charges.) 

Table [I] presents results for linear (A = 0), rotating (afl — 0.3) models, for scalar source points with unit charge. 
For all models the inner boundary conditions were those of the small-x point approximation given in Eq. I|54|l at 
some Xmin- The number of multipoles kept is specified by the parameter ^ max . Choosing £ max = 3 means that modes 
corresponding to monopole and quadrupole were kept; £ max = 5 means that in addition the hexadecapole was kept; 
and so forth. The accuracy criterion used is the quantity "fQ e s/Q, the value of which is unity in the exact solution. 

The results in the table are divided into three sections. In the first section the number of grid points n x , no, and 
n$, was varied, while the values of Xmin and ^ max are kept fixed. The results show 3% accuracy, and demonstrate 
that, for the parameters of this computation there is no advantage to grid size larger than 8001 x 16 x 32. Note that 
simple considerations of truncation error do not apply, since the angular grid is not used in a straightforward finite 
differencing, but rather to establish the angular eigenvectors. 

In the second section the results show that increasing ^ max , for an adequately large grid, improves accuracy, and 
results come within a fraction of a percent of the correct answer. Note that using all the eigenvectors is equivalent 
to no multipole filtering. In that case we would be simply doing finite differencing in the multipole basis, and we 
would be plagued by the problems described at the start of Sec. IIIII Accuracy must, therefore, drop off when £ max is 
increased past some optimal value. The results in the table suggest that the optimal value for this model and this grid 
size may be £ max ~ 5. Larger values of ^ max are more difficult computationally, and appear to give no improvement 
in accuracy. 

In the third section of the table the value of Xmin is decreased and, for a fairly large grid and for £ max = 5 the results 
show that the errors in the inner Dirichlet data were the dominant source of error. More important, the results show 
that very high accuracy can be achieved with the eigenspectral method using a small number of multipoles. 

In the last column of table ^ several results are given of a "two region" method of computation. The motivation 
for this method is that severe multipole truncation is really necessary only in the wave zone. Closer to the source 
more mutlipoles can be kept and more precise computation can be carried out. For the results in the last column, 
all multipoles up to ^ max were kept for grid points with x < 3a; for \ > 3a only the monopole and quadrupole 
eigenmodes were used. The results show no increase in error compared with standard method, but the error in any 
case is dominated by the inner boundary data, not by truncation. 

For nonlinear models, with A = —25 and ^0 = 0.15, table ITU gives results roughly equivalent to the linear-model 
results in table [lj Now there is no a priori correct answer known, so we look only for convergence of the value of 
the monopole moment Q c g. (Due to the effects of the nonlinearity, this value can be reduced well below unity.) 
The computational results in the table show few differences from those in table [lj Again, the answer is shown to be 
stable for moderate grid size, and there is no evidence of a strong dependence on £ max . The two-region computations 
converged more quickly than those with a uniform multipole cutoff, and give results in good agreement with those of 
the uniform cutoff standard approach. This two-region technique, therefore, can be considered a computational tool 
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that may prove useful in more difficult problems. 

Though there is no a priori known general solution for the nonlinear problem, we do know one useful limit of the 
solution. As argued in Sec. IIVI and in Paper I, should have an approximately Yukawa form for a range of small 
X- Evidence of this in the results is presented in Fig. EI which gives computed nonlinear outgoing solutions near the 
sources. The computations start with the boundary conditions of Eq. (|54|l . The variable R in the figure is Z — a along 
a line through the sources, that is, the radial distance from a source. A straight line in the log-log plots of the figure 
indicate that 'J is falling off approximately as 1/AirR; the downard deviation from a straight line is a manifestation 
of the nonlinearity. According to the analysis following Eq. (|49|l . the radius i?ii n at which nonlinear effects become 
significant, decreases with increasing |A| and with decreasing v&q- Results for our standard choice ^>q = 0.15 are shown 
on the left. The nonlinear effects become important for R/a=Ru n /a on the order a few tenths. In this case \& is 
comparable to when nonlinear effects become important, and the Yukawa form is not distinguishable from a l/R 
fall off. For more convincing evidence of the working of the nonlinearity we change Wo to 0.01. The results, shown on 
the right in Fig. [21 for A = —25 shows the excellent agreement of of the computed solution to the Yukawa form in the 
range R/a=0.1 to around 0.4. 




FIG. 9: Near-source fields for nonlinear models. 

Paper I used Eq. (|50|l as the basis of an estimate of the nonlinear solution. That estimate was applied to the radiation 
"reduction factor," the factor by which the radiation amplitude is reduced for a nonlinear model as compared with 
a model with the same parameters, but with A — 0. This provides us with a convenient comparison of the nonlinear 
results in Paper I and with the present eigenspectral method. In table lllll we give those Paper I results again, along with 
eigenspectral computations of the same models. We present additional models, since the standard coordinate/finite 
differencing method of Paper I was limited in the size of A for which Ncwton-Raphson runs converged; with the 
eigenspectral method we can give results for much larger values of —A. The last column in table llllj gives the computed 
reduction factor computed with the eigenspectral method keeping only the monopole and quadrupole terms. (Note: 
In Paper I the factor 1/7=1/1.048 was mistakenly omitted from the delta function source. Here we choose to treat 
that as a source of strength 1.048, rather than unity. Our eigenspectral computations therefore used this enhanced 
source. The estimates of Rn n were also slightly in error in Paper I since they assumed a unit source strength. They 
have been recomputed and are slightly different from the estimates presented in Paper I.) 

The agreement of the computed results with the simple estimate is gratifying, as it was in Paper I. More important, 
the comparison of the second and third columns of tabic IIIII shows that the eigenspectral method with only two 
multipoles gives 1% agreement of the computed radiation with the very different and much more computationally 
intensive finite difference method. 

We have argued that the details of the higher moments of the source are not important in determining the radiation. 
Some numerical justification for this is given in Fig. 1101 which shows computed results for nonlinear models with 
the standard parameters. The solid curve uses the point source initial data of Eq. I|54|) as inner Dirichlet data 
at Xmin = 0.2a. For these inner boundary conditions the multipole moments at Xmin = 0.2a are a = —13.90, 
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TABLE III: The radiation reduction factor due to the nonlinearity. For all cases, = 0.15, and aQ = 0.3 . The second column 
refers to Eq. 1501 . The third column gives the reduction factors presented in Paper I. The last column gives the results of the 
eigenspectral computation with Xmin = 0.3a, with outgoing boundary conditions at Xmax = 50a, i m&K = 3, and a grid with 
n x — 8001, ne = 16, n$ = 32. 



A 


Estimate 


Paper I 


Eigenspec 


-1 


69% 


78% 


77.0% 


-2 


62% 


68% 


68.1% 


-5 


53% 


55% 


55.7% 


-10 


46% 


47% 


46.4% 


-25 


37% 


35% 


35.5% 


-50 


25.6% 




28.5% 


-100 


19.5% 




22.7% 



0.001 




r/ a 



FIG. 10: For outgoing nonlinear waves, the sensitivity of the radiation to details of source multipole structure. 



O20 = 0.11045. The coefficient corresponding to the real part of Y22 is 0.1920; the coefficient corresponding to the 
imaginary part is zero. We first compute the outgoing linear solution for these inner Dirichlet data. Next we, somewhat 
arbitrarily, set all the quadrupole components to -2.209, which is 20 times the original value of 020, and calculate 
the outgoing linear radiation. The results in Fig. ^| show that the effect on the radiation is of order 10%. Some 
interpretation is needed to connect this result to multipoles of sources, especially because the effect on the radiation 
of a physically plausible source quadrupole depends on the size of the source. 

If we had a source with a surface at x = 0.2a the computed result tells us that a rather large deformation, with 
|<22fc/io| ~ 0.16 will have a 10% effect on the radiation as compared with a source with a negligible quadrupole. 
With a simple argument, we can apply this 10% effect to sources of other size. Mathematically the conditions at 
Xmia = 0.2a can be ascribed to a source with a surface at Xsurf 0.2a. Quadrupole moments fall off as 1/-R 3 , 
where R is the distance (small compared to a) from the source point, and monopole moments fall off as 1/R. The 
quadrupole to monopole ratio, therefore, falls off as 1/R 2 , or 1/x 4 - Since the 10% effect corresponds to \a2k/do\ ~ 0.16 
at Xsurf = 0.2a, it also applies, to 0.16 x (1/2) 4 w 0.01 at Xsurf = 0.4a, and to 0.16 x (2) 4 w 2.6 at Xsurf = 0.1a. 
This means that the radiation generated is reasonably sensitive to a mild quadrupole deformation of a source that 
is comparable to the size of the binary system, but for small sources unphysically large deformations are required to 
have any effect on the radiation. (See also the related discussion of the two-dimensional case in Appendix [BJ) 

Figure 1111 shows the central result of our method, the accuracy of the outgoing approximation, for a model with 
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il = 0.3 and A = —25, one of the models presented in Paper I|10j. A measure of the strength of the nonlinearity 
is the fact that the nonlinearity reduces the amplitude of the waves to 35% of those for A = in the same model 
(i.e. , the same fl and inner boundary data and outgoing boundary conditions). The figure shows that the extracted 
solution is in remarkably good agreement with the computed nonlinear solution in the three regions of the extraction 
protocol described in Sec. IIVI (i) the wave region in which the solution is treated as a half-outgoing and half-ingoing 
superposition, (ii) the inner region in which the outgoing solution is taken to be well approximated by the standing 
wave solution since the radiative boundary conditions are irrelevant close to the source, and (iii) the blending region 
described in Eqs. ||S0jH|nD- 

The excellent agreement between the computed outgoing solution and the extracted approximation should not be 
confused with agreement with the exact solution. As we have seen in the comparisons of exact and computed linear 
solutions, e.g. in Fig. |SJ the mutlipole filtering does entail an inaccuracy of a percent or so. Figure ITT1 then, is not a 
demonstration of the accuracy of the eigenspectral method, but rather a powerful statement about effective linearity, 
the accuracy of the process of extracting an outgoing approximation from a standing wave solution. 




FIG. 11: Comparison of a computed outgoing nonlinear solution (continuous curve), and an approximation to the outgoing 
solution extracted from the standing wave solution (data type points). Results are shown for a typical nonlinear scalar model, 
with parameters f2 = 0.3, A = — 25, = 0.15, Xmin/a = 0.05, Xmax/a = 200, for a grid with n x x tie x n$=12001 x 16 x 32. 
Results are shown along the Z axis. The solution is plotted as a function of Z, the distance, from the origin, along the axis 
through the source. The extracted points in the wave zone are a result of treating the waves as linear. The small-distance 
plot shows the blending region and the inner region in which the standing wave solution is used as an approximation for the 
outgoing solution. 

In Table IIVI we present a broad overview of the validity of effective linearity for a range of nonlinear strengths 
in aQ, = 0.3 models. As in table IIIII we give the nonlinear "reduction factor," the reduction in wave amplitude 
due to nonlinear effects. Here we present a comparison of those factors for computed outgoing solutions and for the 
approximate outgoing solution extracted from the standing wave solution. It is clear that judged by the criterion of 
reduction factor (and limited to afi = 0.3 models), effective linearity is highly accurate, within a percent or so, for 
models with extremely strong nonlinear effects. In the = 0.01, A = —100 model, the nonlinearity reduces the wave 
amplitude by a factor of 40, but effective linearity appears to be accurate to better than 1%. It should be noted that 
the agreement of the computed outgoing solution and the extracted outgoing solution is excellent even for models 
(e.g. , "to = 0.01 and small A) for which the strong nonlinear effects are not confined to a small region around the 
source points. This is evidence that effective linearity does not require such confinement; it only requires that the 
nonlinearity falls off before outer boundary effects are important, i.e. , in the induction and wave zones. 



VI. CONCLUSIONS 



The fundamental concepts of the PSW method were introduced in Paper I. In the current paper we concentrate on 
efficient numerical methods for solving the mixed PDEs of the PSW method. The innovative method we present here 
is a mixture of adapted coordinates, multipole filtering, and the use of eigenvectors in place of continuum multipoles. 
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TABLE IV: The reduction factors for nonlinear outgoing waves (the decrease in amplitude due to nonlinear effects). For 
aQ — 0.3 models, the factors are compared for the directly computed outgoing solutions and for outgoing solutions extracted 
from nonlinear standing waves solutions. The A = results indicate linear models in which the reduction factor is unity by 
definition. The value 1.0064 found for the extracted solution gives an indication of the numerical accuracy of the extraction 
procedure. All results were computed with ^ max = 3, Xmin = 0.2a, and Xmax = 50a on a grid with n x x ne x n$ = 8001 x 16 x 32. 
The reduction factor was computed by taking the ratio of the quadrupole components. Also listed are the estimated values of 
i?iin, the distance from the sources beyond which the nonlinear effects are suppressed, and estimates of the reduction factors 
based on the estimates of Ru n - [See Eq. 11501 1 For consistency with table ITTT1 the source strength has been taken to be 1.048. 
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This method seems to meet the needs of the problem remarkably well. The method requires relatively little machine 
memory, and runs very quickly on workstations. The power of the method has allowed us to run nonlinear scalar field 
models with larger velocity and much larger nonlinearity than was possible with the method of Paper I. 

We have shown that the method is convergent and reliable in a number of senses: (i) For a linear problem, the 
computed solution converges to the known analytic solution as the computational grid becomes finer and the number 
of retained multipoles increases, (ii) For a nonlinear model the Newton-Raphson iteration stably and reliably gives a 
solution to outgoing or standing wave problem. We have confirmed that our solutions agree, to the expected accuracy, 
with the results presented in Paper I. 

In addition to the role they play in the efficient computation, the adapted coordinates are very well suited to the 
specifications of inner boundary conditions, rather than to the specification of actual source terms. We have confirmed 
that there is low sensitivity to the details of the inner boundary conditions. The solution in the wave zone has a 
sensitivity to these conditions that is compatible with physical intuition; there is no excess sensitivity that is an 
artifact of the numerical method. 

Two major points are worth emphasizing. First, we have confirmed that excellent results can be obtained for 
moderate source velocities with computations that keep only the monopole and the quadrupole moments of the 
adapted coordinates. This allows an enormous decrease in the computational intensity of a solution. The cost is only 
a moderate increase in analytic complexity. A second and even more important point concerns "effective linearity," 
the approximate validity of superposing half- ingoing and half-outgoing nonlinear solutions. We have been able to 
verify effective linearity for a wider range of nonlinear models than in Paper I, including models with extremely strong 
nonlinearity. 
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APPENDIX A: COEFFICIENTS FOR ADAPTED COORDINATES 

The inner products V% • V0, Vx ■ V<&, and V0 • V$, vanish since the adapted coordinates are orthogonal (with 
respect to a Cartesian metric on X,Y,Z). The other inner products and Laplacians are evaluated with the explicit 
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transformations in Eqs. H19fl - I|24[) . from which we find 
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where Q is the function 



Q = v /a 4 + 2a 2 X 2 cos(29) 



In general, the A and _B terms are computed from the following: 
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In the case of the TCBC coordinates defined in Eqs. I|19l) - (|24|l . the explicit forms of the coefficients are 

- _ a 4 sin 2 (29) cos 2 $ 



(A4) 
(A5) 
(A6) 

(A7) 

(A8) 
(A9) 
(A10) 
(All) 
(A12) 
(A13) 
(A14) 
(A15) 
(A16) 

(A17) 



A 0e 



COS' 



$ [x 



2 1 a 2 cos(29)]' 



(A18) 



= sin 2 $ 



X 2 cos(29) 



Q - a 2 - x 2 cos(29) 



(A19) 



23 



a 2 [x 2 + a 2 cos(28)] sin(28) cos 2 $ 



(A20) 



.4 



[Q + a 2 + x 2 cos(20)] sin $ cos $ 



(A21) 



^4e* — — 



sin($) cos($) [a 2 + X 2 cos(2 8) + Q] [ X 2 + a 2 cos(2 8)] 
X 4 sin(2 8) 



(A22) 



5 X = 



[cos 2 ($) {3a 2 cos 2 (28) - Q - 2a 2 + X 2 cos(28)} + Q + a 2 + X 2 cos(28)] 



(A23) 



(3Q + a 2 + x 2 cos 28) sin($) cos($) 
Q-a 2 -x 2 cos 28 



(A24) 



where 



VQ + a 2 + X 2 co S (28) (ccos2$ + d) 



XVQ"« 2 -X 2 cos(28) 



c = a 2 x 4 cos(2 8) + 2 a 4 X 2 + 4 a 6 cos(2 8) + 4 aV (cos(2 8)) 2 - 4 a 4 Q cos(2 8) - 2 a 2 Qx 2 - X 6 



d = x 4 (a 2 cos(2 8) + x 2 ) ■ 

The coefficients needed in the Sommerfeld boundary condition Eq. (|27|) are 
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APPENDIX B: THE STANDARD SPECTRAL METHOD FOR THE 2+1 DIMENSIONAL LINEAR 

SCALAR FIELD 



Here we consider the 2+1 dimensional version of our helical problem, one equivalent to the 3+1 problem with line 
sources that are infinitely long in the z (equivalently Y ) direction. We choose to set A = 0, i.e. , to make the problem 
linear, since that will turn out to allow a very efficient method of multipole projection. The 2+1 dimensional version 
of Eqs. © and (JSU is 
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Here ip = <f> — fit where <f> is the usual polar angle tan 1 (y/x) in the x, y plane. 
The 2+1 dimensional forms of the adapted coordinates of Sec. ITTlare 
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With these coordinates the 2+1 dimensional version of Eq. takes a form like that of Eq. (jSJ. As in the 3-dimensional 
case we use only the homogeneous form of Eq. HB1|) . since the effect of the source will be introduced through inner 
boundary conditions. 

In working with the sourceless linear 2+1 dimensional problem, it turns out to be convenient to divide the original 
wave equation by Q/y 2 , where 
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By dividing through by Q we have put the wave equation in a form in which the coefficients are 9-independent in the 
f7 — ^ limit. With the standard method, used below, for projecting out multipole components of the wave equation, 
this property of the coefficients means that the mixing of the multipoles can be directly ascribed to the rotation. 
We now expand the standing wave solution ^(x, 9) as 
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and our equation becomes 
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Projecting with J Q * sinm9 • • • dQ, gives zero by symmetry; projecting with J Q * cos m9 • • • dQ, gives 
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When the explicit expressions for A-E, in Eqs. HB6(1 ~ I|B10(1 . are used in Eqs. (|B14fl — (|B16|) . the results are 

4 f n / 2 sin 2 (29) 
m„,„ (,,,,, - ti 2 a 4 — I cosm9cosn9 — dQ 



(B15) 



(B16) 



(B17) 



_n2 4^2 

X 7T 



,| I/2 n Q (a 2 cos(29) + x 2 ) 2 
n / cosm9cosn9 — a9 



r /2 . 2a 2 (2a 2 cos(2 9) +x 2 )sin(2 9) 

— n I cosm9sinn9 

/o 



Q 



(B18) 



_ 1 4 a 2 fl 2 



ir/2 



cosm9 cosrt9- 



+ 3 a 2 cos 2 (2 9) + 2 x 2 cos (2 9)) 

Q 



de 



r/2 _ S in(29) (a 2 cos(29)+ X 2 ) 

-2n / cosm9smrt9 — — a9 

Jo Q 



(B19) 



where 



2 if to = n = 

< „,.„ = 5 m +n.O + #m,n = { 1 if TO = TV ^ 

if to 7^ n 



(B20) 



The integrals needed in Eqs. (|B17|) (|B19|I are all of the form 
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where P and J are integers. With trigonometric identities, and with the substitution x — sin 9, such integrals can all 
be expressed in terms of integrals of the form 



JC N (k) 



,.N 



where 



o vl — fc 2 ^ 2 -^! — x 2 
2a X 



dx , 



k = 



a 2 + x 2 



(B22) 



(B23) 



For the integrals in Eq. i|B21(l only even values of N are needed for /Cat. All such integrals can be evaluated in terms 
of the complete elliptic integrals |l8|, 
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To use the elliptic integrals to evaluate the fCj^(k), for even TV, we start with the relationship (for M > 0) 
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This gives us the recursion relation 

(2M + l)/C 2M (fc) - (2M + 2)(1 + k 2 )JC 2M +2{k) + (2M + 3)fc 2 /C 2A f +4 (fc) = 

We know that 



(B25) 
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K {k)=K{k) , 



and we can easily show that 



k 2 ih = ±[K(k)-E(k)} 



so all values of /C2m(&) follow from the recursion relation. 

Very efficient computation follows from the results above. At a given value of x, only two integral evaluations must 
be done, those for K(k) and E(k). All values of /C 2 M(fc) then follow from the recursion relation, and hence all values 
of a mn , Pmn, Iran-, can be found at negligible computational expense. 

Here we use this method only to illustrate important general issues. To make this illustration as clear as possible, 
we take the simplest nontrivial case of the multipole expansion: we keep only the n — and n = 2 terms [T^. so that 



*(x,e) = a (x)+a 2 (x) cos 26 



(B27) 



The equations in Eq. I|B13[1 . then describe the interaction of the monopole ao(x) and quadrupole a 2 (x) terms. For 
x/a -C 1 and x/ a 3> 1 the terms that mix the multipoles die off; it is only in the transaction region, \j a ~ 1 that 
there is strong mixing of the multipoles, a mixing that for our problem is quadratic in the source velocity aft. The 
process of the generation of radiation can be viewed as the growth of a 2 from the small- x near-source region to the 
large- x radiation region. 

In principle, the linear standing wave problem could be solved by including a term a s sin(26), leading to an 
additional second order differential equation. At some inner boundary Xmin a the values could be specified for all 
multipoles, and at some outer boundary Xmin 3> a, a fall off condition could be specified for do- Ingoing or outgoing 
conditions could be used to relate a 2 and a s . 

It is easier, and more instructive, to use another approach to finding the standing wave solution, the minimum 
amplitude method presented in Paper I. For the linear problem it is straightforward to show[j| that of all solutions 
that (i) have the form of Eq. I|B11|I . and (ii) correctly couple to the source, the solution with the minimum wave 
amplitude in each multipole is the standing wave solution, i.e. , the solution that is half-ingoing and half-outgoing. 
In principle, the minimum-amplitude criterion can be used as a definition of standing waves in a nonlinear problem. 
Here we are dealing with a linear problem, so we are simply exploiting a known property of the solutions. 

In this minimum-amplitude method we specify ao, dao/dx, a 2 and da 2 /dx at Xmim then shoot outward. The choice 
of do and dao/dx, at Xmin are those for unit point charges. The value of ao at Xmin sets the scale of the linear solution. 
The value of dao/dx can be approximated as l/irXnua- I n principle this starting value can be adjusted so that at 
large x the results for ao satisfy the fall-off condition dag/dx = ao/xl°SX- I n practice, the overall result (after the 
minimization described below) is very insensitive to the choice of dao/dx at Xmin- 

The choices for ao and dao/dx at Xmin are determined by the form of \P very close to one of the unit point charges. 
Some care is necessary in taking this limit. In the x,y frame comoving with one of the unit point charges the value 
of $ due only to that charge is log \(x 2 + y 2 ) / a 2 ] / 4ir . We must now transform this result to the "lab" frame of our 
calculation in which we use coordinates x, y, t. The Lorentz transformation is 



y = iv- 



(B28) 



(The last relation follows since y is in a coordinate frame that comoves with the point source, but which is related 
to the lab frame by a simple translation by vt.) Since $ is a Lorentz invariant we have \& = log (r 2 /a 2 )/47r where 



= x + 7 y , expressed in adopted coordinates, is 



It. 

4 a 2 



lX 2 
2 a 



cos (26) 



+ ( 7 2 -l) \ sin 2 (26) \ 
'A a A 



\ cos 26 
a 2 



(B29) 
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and \& is therefore 



^ = 4^ * l0g 



4a 4 



log[l + ( 7 2 -l) sin 2 (26)] +0(x/a) 



l0g 



4a 4 



2 log 



7 + 1 



0(cos (49)) + 0( X /af 



(B30) 



The last relationship is meant to emphasize that in the \ limit 4" has no cos (29) component. This comes from 
the fact that the Lorentzian pancaking of the source field has the nature of a local quadrupole deformation, while the 
cos (29) term represents a local dipole. (As defined in Eq. (|B3(1 . 9 « 8/2 or (8 + it)/ 2 near the source points, so the 
cos (29) dependence near the source point corresponds to the dipole field of the source.) 

The contribution, near x = 0, due to the distant point can be approximated with r — 2a so that at x — Xmin the 
appropriate starting conditions for unit point charges are 



X 

4a 4 



2 log 



7+1 



, da 1 
log4 — = — 

dx kx 



n da2 n 
a 2 = — — = 

dx 



(B31) 



While dai/dx at x ~> vanishes in principle, in the numerical computation, dai/dx plays a more delicate role. It is 
chosen to minimize the wave amplitude at large x- It is actually the values of a 2 and dai/dx, at Xmin that determine 
the radiation field at large Xi and determine whether there is any a 2 radiation except the radiation coupled to the 
source. To minimize that radiation (and suppress radiation that is not coupled to the source) we fix a 2 and vary 
da-i/dx- (We could just as well fix da 2 /dx and vary a 2 .) Minimization is taken to mean the minimum value of the 
amplitude defined by 



(B32) 




Amps (x/a)VMo» + + -± 



For the expected large distance form of the waves r -1 / 2 cos (2flr + 8), the quantity inside the angle brackets is expected 
to be nearly ^-independent; this is confirmed by the numerical results. The angle brackets denote an average over a 
wavelength, resulting in a quantity that is x independent to high accuracy at large x- In our minimization procedure 
we vary da^/dx at Xmin to minimize Amp. The meaningfulness of this minimization procedure, of course, depends 
on its insensitivity to the details of how a 2 and c?a 2 /dx sxe ch.0S6n 3,t Xmin- 

This is discussed below. 

1.5 i . , . , . , . 1 
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FIG. 12: Comparison of the analytic waves and the computed waves for afi = 0.5, Xmin = 0.05a and Xmax = 70a, and ai = 
at Xmin = 0. The monopole is subtracted in the figure on the right to allow for a comparison of computed and analytic wave 
amplitudes. . 



For the source speed a£l = 0.5, results, for 9 = 0, of our computation are shown in Fig. ^1 and are compared with 
the analytic solution for two point sources each of unit strength 



* = Vl -a 2 fl 2 



— In (r/a) for r > a 
for r < a 



J m (mr2r<)A r m (mSlr>) cos (m9) 



m=2,4,6... 



(B33) 
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TABLE V: The errors (differences from analytic solution) of the monopole+quadrupole approximation for different values of 
the source speed aQ. The overall error is the difference of the computational and analytic solution at Xmax = 70a; the wave 
error is the difference of the computational and analytic peak-to-peak amplitudes. 





Overall Error 


Wave Error 


0.1 


0.08% 


0.14% 


0.3 


1.0% 


0.1% 


0.5 


3.9% 


1.2% 


0.7 


13% 


3% 



(where r>, r< indicate, respectively, the greater and lesser of r, a). For our more typical choice ail = 0.3, the difference 
of the computed and analytic solution are too small to show up well in a plot. The plot on the left shows the comparison 
for the whole range of the computation. Dividing the analytic solution by 1.039 brings it into nearly perfect alignment 
with the computed solution; we therefore characterize the overall error in the computed solution as 3.9%. Overall 
errors for several values of ail are listed in Table |V| 

The plot on the right in Fig. El focuses on the oscillations in the far zone by removing effective monopole terms. 
For the analytic result this is done by subtracting the asymptotic monopole solution VT — a?il 2 log (x 2 /2a 2 )/7r. For 
the computed solution this is done by plotting only a2(x). The comparison shows two effects. First, there is a 1.2% 
difference in peak-to-peak amplitude. Errors of this type are tabulated for different values of ail in Table A 
second and more interesting effect, apparent in the waves of Fig. 1121 is the difference in shapes. Since the computed 
wave contains only the quadrupole component it has a nearly perfect sinusoidal form. The analytic solution, on the 
other hand, shows a rapid rise and a slow fall off due to the contribution of the hexadecapole and higher modes. 
For small ail the amplitude of a cos(m9) component depends on ail approximately according to (ail) m ~ ' 5 , so the 
contributions from m = 4, 6, 8 . . . decrease quickly with il. 

We now want to use the accuracy of the linear 2+1 dimensional model to make an important point about the 
role played by the source structure in generating radiation. To do this in our monopole+quadrupole computational 
models we depart from considering "point" sources and we vary the choice of a2(xmin). In each case the specification 
of daijdx at (xmin) is fine-tuned to get a minimum wave amplitude at Xmax- For this investigation to have physical 
relevance we need to ask what a "reasonable" value is for (^(Xmin)- 

If Xsurf specified the (approximately spherical) location of the outer surface of the source, then we can write 

02(Xsurf) = Ka (Xsurf) , (B34) 

where k indicates the relative quadrupole strength of the source. We expect k to be small for realistic sources and of 
order unity only for highly distorted sources. 

In our computational models we take Xmin/a to be very small, typically 0.05. It is useful, however, to consider the 
computed results applying to larger sources, sources with Xsurf significantly larger than Xmin- To do this we can use 
the fact that outside, but very close to a source, ao varies as logx and a 2 falls off as 1/x 2 - From this we conclude 
that for small x 

a 2(Xmin) = k Xsurf fog Xsurf 
ao(Xmin) Xmin l0 g 

Xmin 

We can now fix K = ±1, its maximum reasonable value, and we can choose a value of Xsurf, the value at which 
the quadrupole and monopole have equal strength. These choices determine a2(Xmin)/ao(Xmin), and therefore the 
computational model. 

Results are presented in table IVTl For the choices Xsurf/a = 0.1, 0.2 and 0.3, and for k =1 and -1, the amplitude 
of the quadrupole waves is computed for Xmin/a = 0.05 and Xmax/a = 70. For each choice of a 2 (xmin), the starting 
value of daijdx is fine tuned for minimum wave amplitude. Since a2(x) should have the form const. /x 2 , the value of 
daijdx should be very nearly equal to —2ai/x- The values for daijdx in table IVll are very close to this prediction. 
The table presents the "relative amplitude," the ratio of the radiation amplitude to the radiation amplitude for the 
case a2(xmin) = 0. It is clear from the results that the structure of the source has little influence on the radiation 
unless the source is large and has an extreme nonsphcrical deformation. 

In the case of gravitational sources, an even stronger statement can be made. As already point out, the "quadrupole" 
mode a 2 (x) cos(20) near the source point is actually a local dipole ai (x) cos(#i). This dipole moment can be viewed 
as a displacement of the center of scalar charge, and hence a change in the radius at which the source points move. 
The change in the radiation amplitude can be ascribed to this change in radius. This explanation of the change can 
be made quantitative. We let r represent the distance from source point 1 at which inner boundary data is being 
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TABLE VI: The effect on the wave amplitude of the conditions on a?, and da-zjdx at Xmin = 0.3a. The role of the ai term can 
be understood as an effective shift of the center of scalar charge. See text for details. 



K 


Xsurf 


ffl2(Xmin) 


dai 1 d\ at Xmin 




1 — 4z7VVO,2 f & 


l 


0.1 


-2.92007 


116.999652 


1.045 


1.046 


-l 


0.1 


2.92007 


-117.011864 


0.955 


0.954 


l 


0.2 


-8.16416 


327.127955 


1.125 


1.128 


-l 


0.2 


8.16416 


-327.140167 


0.875 


0.872 


l 


0.3 


-13.7416 


550.612310 


1.211 


1.216 


l 


0.3 


13.7416 


-550.6245521 


0.789 


0.784 



specified and we let 5 be the radial distance by which the center of scalar charge is being moved outward. The solution 
for the scalar field is 

* = —log (r 2 - 2r6cos6 1 + S 2 ) « — [log r 2 - 2(5 /r) cos 6 A . (B36) 
in v ' 47r L J 

From this we infer that 

02=-/-. (B37) 

This means that the radius of orbital motion is changed from a according to 

a — > a + 5 = a — 2-Krai . (B38) 

From Eq. (|B33|1 we see that for quadrupole radiation (that is, m = 2) the amplitude of the waves scales a a 2 . It 
follows that the amplitude of the waves should depend on a 2 according to 

amplitude oc 1 — inra-ij a . (B39) 

For the computations presented in table I VII a 2 is evaluated at Xmin = 0.05a, so r ~ 0.05 2 a. The numerical results 
following from this simple explanation of the shift of the center of charge, presented as the last column in tabic IVll 
are convincingly accurate. 

This explanation for the role of the m = 2 mode is of considerable significance for gravitational problems. The 
equivalence principle implies that there is no local dipole for the gravitational sources. Thus the starting value of a 2 
which, at large distances gives the quadrupole radiation, is not a parameter of the structure of the source; if we know 
the location of the effective center of the source, the structure is fixed. 



APPENDIX C: DETAILS OF THE EIGENSPECTRAL METHOD 



In this appendix we explain how the continuum angular Laplacian of Eq. I)35|l is implemented as a linear operator 
in the N — n@ x n$ dimensional space. That linear operator must represent the angular Laplacian evaluated at a 
grid point o , That is, the linear operator L a &,r/ must satisfy 

[sineV^*]^«£La6,y*«, (CI) 

where the approximation is due to FDM truncation error. It will be convenient below to write the linear operator as 
a sum L = IS 1 ' + with containing 9 derivatives and L^> containing $ derivatives. 

By exploiting the symmetries of the PSW configuration we can limit the range of angular coordinates to one 
quarter of the complete 2-sphere. The indices a, i range from 1 to ne , representing, respectively Q — A0/2 to 
O = 7r/2 — AO/2. The indices b, j range from 1 to n$, representing, respectively <& = A<I>/2 to <& = it — A$/2. Our 
goal here will be to show that with this choice of the grid, the linear operator can be chosen to have the symmetry 

L a b,ij — Lij ja b . (C2) 

The elements of L^ ij have a different form for the boundaries at O = AO/2, ir/2 — AO/2 and at <I> = A<f>/2, 
7r — A$/2, and in the interior of the angular grid. We consider each case separately. 
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Case I: Interior points, 1 < a < ne, 1 < b < nj> 



The first contribution is 



L ab,ij = T^oy [ sin0 a-i/2^a-i,i - (sin Q a +i/2 + sin G Q _ i /2 ) S ia +sin8 a+ i/ 2 ^o+i,i] (C3) 



and the second is 



Here (5 JjC is the Kronecker delta, and sinG a ±i/2 is defined to mean sin (Q a ± AS/2). It can be seen that both 
contributions to L a b.ij are symmetric with respect to the interchange of the pair ab with the pair ij, and hence 
Eq. ljC2() is satisfied. 



Case II: Boundary at a = 1, 1 < b < na> 



The derivative part of the operator formally takes the form 

(1) _ 5 bj 



L 



lb,ij 



(AG) 2 



[sin 8o t i — (sin + sin A6) 5n + sin A6<5 2 , 



(C5) 



In the sum in Eq. I|C1|I the So,i term respresents ^(G = — AG/2, $), which is not a value available on the angular 
grid. This term however, is multiplied by sin = and can be ignored, so satisfies the symmetry condition in 

Eq. l)C2Jl . The form of ^ also applies without change to a — 1, and hence Eq. I)C2J| is satisfied for the index range 
a = 1, 1 < b < n$. 
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Tl+AO/2 f 

jc-AO/2 r - - - 
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FIG. 13: An angular grid with ne = 5 and n$ = 6. Grid points, points at which a value for ^ is computed, are connected by 
solid lines. The dashed lines extend the grid to "phantom" points needed for the computation. For the FDM implementation 
of the Laplacian at point A the value of "J at point A' is needed. By the symmetry of the physical problem, this value can 
be replaced by the value at point B, which is on the grid. Similarly the value at B' , when needed, can be replaced by that at 
point A; the value at C' can be replaced by that at C; the value at D' can be replaced by that at D; and so forth. 



Case III: Boundary at a — ne, 1 < 6 < n$ 
In this case the G derivative part of the operator formally takes the form 

L nlh,H = 7^1)2 t sin - - ( sin (*/2 ~ AG) + sin (tt/2)) 5 ne4 + sin (n/2)6 ng+hl ] . (C6) 
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The 8 ns +\ t i represents fy(Q ne +i, = \&(0 = tt/2 + A0/2,$). This value is not directly available on the grid, but 
we can get an equivalent value that is on the grid by using the symmetry 

#(e,$) = *(7r-e,7r-$). (C7) 

This is equivalent to the statement that VE 1 is invariant with respect to inversion through the origin and reflection in 
the orbital plane (or equivalently Z — ► —Z, X — > —X, Y — » Z in Fig. [3J) As shown in Fig. we can, therefore, 
replace *(6 = tt/2 + A6/2,$) with *(6 = tt/2 - A0/2, tt- $), or can replace f(9 ne+1 , $ 6 ) with *(0„ e ,^) where 
6 = A$ — 6. Equivalently, we can rewrite Eq 1C6(1 as 



(A0) 2 L— v/- r-r,- / ■ — v/-yy-T.e,»j ■ ^ A qj 

The term that has been introduced is 



L nlb,r 3 = TaHW [ sin (^/ 2 - A©M„ 9 - M - ( sin (tt/2 - A0) + sin (tt/2)) * Beii ] + sin (tt/2)^ . (C8) 



L neb,n e b ( AflV ' ( C9 ) 



(A0) 
symn 

the same as in Case I, and hence Eq. (|C2() is satisfied for the index range a = ne, 1 < 6 < n<j, 



Since 6 = 6, we have L*' 1 '', 7 = L^'r , which satisfies the symmetry in Eq. 1C2II . All other terms in L^K . remain 

' neb, "eb neb, neb * » 19»,'3 



Case TV: Boundaries at 1 < a < ne, & = 1 and 6 = n$ 
For these boundary points Case I considerations apply to L^-. For 6=1, however, takes the form 

The (5o,j refers to an angular location (<& = — A$/2) that is not on the grid. Here we can use the symmetry 
*(0, -$)=*(©,$), and hence *(0, -A$/2)=*(6, A$/2), to replace S 0d with 5\j. The resulting L ab:ij satisfies 
the symmetry of Eq. I|C2|) . 

For b = n$, the considerations are very similar. The \l/(0, — < i>)=5 , (0, <&) symmetry is used to replace 5 niS ,+i.j by 



Case V: Boundaries at a — ne ,6=1 and b = n<j 

Here the forms of - and L„„ - ■ are taken from Case III with 6 = 1, n$, and the forms of l! 2 \ and i^L 
are taken from Case IV with a = ne ■ From the considerations of Case III and Case IV it follows that the results here 
also satisfy Eq. I|C2|I . 

To clarify the results derived above, we list here all nonzero elements of L^ t - and L^^: 



r(i) sine a+ i /2 + sin9 a _i /2 



-T (1) - Sme a-l/2 r I 

^o6,(o-l)6 _ L (a-l)Mi ~~ (A0) 2 



- sin9a + 1 / 2 fn , n<rrio fpi^ 

^ab,(a+l)b - ^(a+l)b,ab - (A6) 2 



L (1) , ,=i (1) , , = 7-4-7 where & = n # -& (C14) 

neb, neb ne6,iie» (AO) 



L Sa> = - sin e a (A9) 2 alll<6< ^ < C15 > 
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1 



(C16) 



sine a (A6) 2 



1 



all b < n$ 



(C17) 



ab,a(b+l) 



a(b+l),ab 



sin9 Q (Ae) 2 



(2) 

a(b— l),afe 



l 



all 1 < b . 



(C18) 



ab,a(b— 1) 



= L 



sin6 a (Ae) 2 



This completes the proof that for the full range of its indices I/ a b,ij satisfies Eq. (IC2|) . With this result in hand we 
can go on to the computation that is central to our eigenspectral method: finding the eigenvectors of 




where the k index indicates that the solution is the kth eigensolution. Aside from the sinO a factor on the right, 
this is a standard eigenproblem for a symmetric real matrix, and we conclude that the eigenvalues are real and the 
eigenvectors form a complete basis. It is easy to show that the factors of sin0 a do not change these conclusions. 

The finite difference problem in Eq. IjCl 9(1 . along with Eqs. IjClj) . can be seen to be the finite difference equivalent 
of the continuum eigenproblem 



With the usual boundary conditions, the solutions of Eq. (|C20(1 can be taken to be the spherical harmonics, and A to 
have values £{£ + 1) where £ is an integer. The solutions of Eq. IjCf 9|) should then be approximately proportional to 
the real and imaginary parts of Yi m (Qi, $j), the approximation becoming perfect as the grid goes to the continuum 
limit. 

We next define the inner product in the grid vector space by the expression in Eq. (|36[) . It is simple to show, following 

(k) 

the usual pattern with eigenproblems, that with respect to this inner product, two nondegenerate eigenvectors Y^ . 

and Y^ are orthogonal as a consequence of the symmetry in Eq. (|C2|) . Since we find the grid multipoles to have no 
degeneracies it follows that the solutions to Eq. I|C19|) constitute a complete, orthogonal basis, and can be normalized 
to satisfy Eq. I|37|) . It should be clear that this is the finite difference equivalent of well known continuum relations. 
In the continuum limit, Eq. (|36|l is the inner product on the two sphere. The orthogonality of our grid multipoles is 
therefore just the finite difference form of the orthogonality of spherical harmonics. 
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v^ g F(e,$) = -Ay(e,$). 



(C20) 



